#ifndef AEROSOL_ATMOSPHERE_H
#define AEROSOL_ATMOSPHERE_H
#include "atmoslayer.h"
#include <QObject>

//parameters used in aerosol distribution caculate
#define R_STEP 0.01
#define R_MIN 0.01
#define R_MAX 10
#define NUM_R_STEPS 1000
#define R_REFREC_INDEX_SEPERATE 1

#define F_AVERAGE_REGION (2*ANGLE_STEP)

//small particals(Aikten partical), 0.01-4 micron
#define R_MIN_1 0.01
#define N_R_1 399

//big particals, 1-100 micron
#define R_MIN_2 1
#define N_R_2 9999

class AerosolAtmosphere : public AtmosLayer
{
    Q_OBJECT
public:
    explicit AerosolAtmosphere(QString layerName,
			       double hBottom,
			       double hTop,
			       double (*smatrixAero)[3][3],
			       double (*pbsdfAeroAveraged)[N_ZENITH][2*N_PI][3][3],
			       double aerosolTau = 0,
			       QObject *parent = 0);
    ~AerosolAtmosphere();

signals:

public slots:


/******************* Aerosol properties**************/
protected:
    double m_aeroTau;
    double m_moleRatio, m_aeroRatio;
    double (*m_smatrixAero)[3][3];
    double (*m_smatrixAeroAndMole)[3][3]; //pbsdf at N scatter angle [0, pi]
    double (*m_averagedFowardPbsdf)[N_ZENITH][2*N_PI][3][3];
    void scatMatrixCaculate();
public:
//caculate aerosol optical depth, scat-Matrix, and save to buffer (*m_smatrix1)[3][3]
    static void aerosolInitTwoPeak(double rMean1,
				   double sigma1,
				   struct complex m1,
				   double rMean2,
				   double sigma2,
				   struct complex m2,
				   double numRatio1_2,
				   double (*smatrix)[3][3]);
    static void readAerosolScatMatrix(QString fileName,
				      double (*smatrix)[3][3]);
    static void averageAerosolMatrix(double (*smatrix)[3][3], double (*pbsdfAero)[N_ZENITH][2*N_PI][3][3]);

protected:
    virtual void solarScatCoeff(int iThetaS, int iPhiS,
				linearPolar *lightOut,
				int scatMode);
    virtual void scatCoeff(int iThetaS, int iPhiS,
			   int iThetaI, int iPhiI,
			   double lightIncident[3],
			   linearPolar *lightOut,
			   int scatMode);

protected:
    //                        zenIn     zenOut    aziOut
    static int (*m_scatAnglesIndexSameDir)[N_ZENITH][2*N_PI];
    static int (*m_scatAnglesIndexDiffDir)[N_ZENITH][2*N_PI];
    static void scatAngleIndexInit();

protected:
    //Stocks rotate transform matrix L1 and L2, each has 2 independent value:
    //L11, L12 for matrix L1; and L21, L22 for matrix L2.
    //symetry L1(thetaI, thetaS, phiI,phiS) = L1(-thetaI, -thetaS, phiI,phiS)
    //same for L2
    //                zenIn     zenOut    aziOut L1  L2
    static double (*m_LSameDir)[N_ZENITH][2*N_PI][2][2];
    static double (*m_LDiffDir)[N_ZENITH][2*N_PI][2][2];
    void LMatrixInit();
};

void saveScatMatrix(QString fileName, double (*pbsdfAero)[3][3]);
#endif // AEROSOL_ATMOSPHERE_H
